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SparseLib++ 


1 About  SparseLibH — h 


SparseLib4-+  is  a C++  class  library  for  efficient  sparse  matrix  computations  across  various  computa- 
tional platforms.  The  software  package  consists  of  matrix  objects  representing  several  sparse  storage 
formats  currently  in  use  (in  this  release:  compressed  row,  compressed  column  and  coordinate  formats), 
providing  basic  functionality  for  managing  sparse  matrices,  together  with  efficient  kernel  mathemati- 
cal operations  (e.g.  sparse  matrix- vector  multiply).  The  Sparse  BLAS  Toolkit  [1]  is  used  to  enhance 
portability  and  performance  across  a wide  range  of  computer  architectures.  Included  in  the  package  are 
various  preconditioners  commonly  used  in  iterative  solvers  for  linear  systems  of  equations.  The  focus 
here  is  on  computational  support  for  iterative  methods  (see  IML++[5]),  but  the  sparse  matrix  objects 
presented  here  can  be  used  in  their  own  right. 

The  goal  of  SparseLib++  is  to  provide  the  ability  to  develop  and  experiment  with  numerical  linear 
algebra  algorithms  through  the  separation  of  the  internal  details  of  a sparse  matrix  representation  from 
the  code  that  uses  it.  With  this  object-oriented  approach,  the  need  for  separate  hand-coded  numerical 
linear  algebra  routines  for  each  sparse  matrix  type  is  reduced. 

The  SparseLib++  library  provides 


• double-precision  sparse  matrix  classes  for  compressed  column,  compressed  row,  and  coordi- 
nate storage  formats 

• access  to  matrix  elements  of  any  sparse  matrix  type  with  the  conventional  notation  A(i,  j) 

• conversion  between  any  two  sparse  matrix  types  through  simple  matrix  assignment  (e.g.  A = B) 

• computational  kernels  based  on  the  Sparse  BLAS  Toolkit  interface  [1]  for  maximum  efficiency  and 
portability  across  various  hardware  platforms 

• basic  preconditioners  useful  in  iterative  methods:  incomplete  LU  (ILU),  incomplete  cholesky  (ICP), 
and  diagonal  scaling 

• the  ability  to  read  and  write  Harwell-Boeing  formatted  files  and  simple  text  files  for  matrix  in- 
put/output 

• easy  integration  with  generic  Fortran,  C arrays  together  with  user-defined  C++  matrix  packages 

2 Sparse  Matrix  Representations 


In  the  following  subsections  we  describe  the  underlying  data  structures  for  the  sparse  matrix  classes 
currently  supported  in  SparseLib++.  Although  the  object-oriented  paradigm  encourages  hiding  such 
issues,  we  have  included  this  section  to  assist  the  user  in  integrating  SparseLib++  with  Fortran  libraries 
and  in  matching  a storage  scheme  with  their  specific  problem  structure.  The  SparseLib++  classes  use  a 
0-based  indexing  internal  representation  (as  opposed  to  the  1-based  indexing  representations  commonly 
associated  with  Fortran,  see  section  5.1).  This  should  be  kept  in  mind  while  reading  the  following 
subsections. 
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SparseLibH — h classes  also  use  the  MVH — [-  matrix/ vector  classes  as  internal  building  blocks.  See  the 
MV++  Reference  Manual[6]  for  details. 

The  current  version  (1-5)  supports  only  double-precision  matrices,  although  it  is  fairly  straightforward 
to  generate  complex  or  other  user-defined  data  type  sparse  matrices  from  the  SparseLib-l  h source  code. 
(We  provide  scripts  to  generate  these.)  Future  releases  of  SparseLib-l — h will  incorporate  fully  templated 
type  parameters  as  template  facilities  mature  in  production  level  C-| — h compilers. 


2.1  Data  Storage  Formats 


To  illustrate  the  various  storage  formats. 

we 

will 

use 

the 

nonsymmetric  matrix: 

/ 
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For  symmetric  matrices,  one  can  store  only  the  upper  (or  lower)  triangular  portion  of  the  matrix.  The 
savings  in  storage  is  traded  for  additional  bookkeeping.  The  current  release  of  SparseLib-j — h does  not 
support  symmetric  storage,  but  future  releases  will  incorporate  this  space-economic  feature. 


2.2  Coordinate  Storage 

The  most  straightforward  scheme  to  denote  a sparse  matrix  simply  records  each  nonzero  entry  together 
with  its  row  and  column  index.  Three  data  structures  are  used:  a val  ( ) array  to  hold  the  floating  point 
values,  and  row  and  column  index  arrays,  rowJ.nd()  and  col_ind().  The  convention  is  that  for  each 
k G {l,...,nz}.  The  value  val (k)  occurs  at  position  (row J.nd(k)  , col_ind(k)). 

The  matrix  A in  (1)  could  thus  be  stored  with  the  following  arrays: 


valO 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

row_indO 

0 

0 

0 

1 

1 

1 

2 

2 

2 

3 

4 

4 

col_Lnd() 

0 

1 

4 

0 

1 

2 

1 

2 

4 

3 

0 

4 

One  should  note  that  the  ordering  of  matrix  elements  in  this  format  is  not  fixed,  and  so  the  given  storage 
representation  for  the  matrix  is  not  unique.  See  section  5.2  for  details. 


2.3  Compressed  Row  Storage 

The  compressed  row  storage  format  views  nonzero  elements  in  each  row  as  a sparse  vector,  storing 
pointers  to  the  first  element  in  each  row  in  row_ptr(),  and  nonzero  values  and  their  associated  column 
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indices  in  the  arrays  val()  and  col^ndO.  An  additional  element  is  appended  to  the  row_ptr()  array 
specifying  the  number  of  nonzero  array  elements. 

The  matrix  A in  (1)  can  be  stored  in  compressed  row  format  with  the  following  arrays: 


row-ptrO 

0 

3 

6 

9 

10 

12 

val() 

1 

2 

3 

4 

5 

6 

7 8 

9 

10 

11 

12 

col_indO 

0 

1 

4 

0 

1 

2 

1 2 

4 

3 

0 

4 

2.4  Compressed  Column  Storage 

The  compressed  column  storage  format  parallels  compressed  row  storage,  but  with  the  roles  of  rows  and 
columns  reversed.  The  appropriate  array  entries  for  the  matrix  A in  (1)  are  as  follows. 


col_ptrO 

0 

3 

6 

8 

9 

12 

valO 

1 

4 

11 

2 

5 

7 

6 

8 

10 

3 

9 

12 

col_indO 

0 

1 

4 

0 

1 

2 

1 

2 

3 

0 

2 

4 

2.5  Creating  Sparse  Matrices  from  C Arrays 

SparseLib++  sparse  matrices  can  easily  be  constructed  from  individual  C/C+4-  vectors.  For  example, 
the  5x5  matrix  A in  (1)  could  be  created  by: 


double  val[12]  = {1 . ,2 . ,3 . ,4 . ,5 . ,6 . ,7 . ,8 . ,9 . , 10 . , 11 . , 12 . }; 
int  colind[12]  ={0,  1,  4,  0,  1,  2,  1,  2,  4,  3,  0,  4}; 
int  rowptr[6]  = {0,  3,  6,  9,  10,  12 

CompRow_Mat_double  R(5 , 5 , 12 , val .rowptr , colind) ; 


See  the  class  man  pages  in  section  A for  details  of  the  various  constructor  parameter  lists  for  each  sparse 
matrix  class. 


3 Sparse  Matrix  Operations 

3.1  Matrix- Vector  Multiplication  via  Sparse  BLAS 


One  of  the  core  operations  in  numerical  linear  algebra  is  the  matrix-vector  multiply,  which,  through 
operator  overloading,  can  be  called  in  SparseLib-f--f  with  the  simple  A*x  notation.  Since  our  goal 
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1 

1 

1 . OOOOOOOOOOOOOOOOe+000 

3 

1 

1 . 6065109611192374e+000 

10 

1 

1 . 1329999260086812e+000 

11 

1 

5 . 54171 56097831336e-001 

13 

1 

1 .49994248007728646-001 

2 

2 

1 .OOOOOOOOOOOOOOOOe+000 

3 

3 

1 .OOOOOOOOOOOOOOOOe+000 

13 

3 

-7 . 0499387778693664e-002 

40 

3 

-4 . 05138316773604806-001 

Figure  1:  An  example  of  a sparse  matrix  text  format:  <i>,  <j>,  <val>. 


in  creating  this  class  library  was  to  allow  not  only  readable  and  reusable,  but  also  efficient  numerical 
linear  algebra  routines,  we  have  provided  matrix- vector  multiplication  operations  which  use  the  proposed 
Fortran  Sparse  BLAS  Toolkit  [1].  The  Toolkit  provides  a uniform  interface  to  kernel  linear  algebra 
routines.  With  vendor  support,  using  this  interface  should  allow  for  efficient  implementations  across 
platforms  by  simply  linking  with  appropriate  machine-specific  BLAS  libraries.  The  current  SparseLib-|— f 
implementation  provides  its  own  sparse  BLAS  routines  which  may  be  used  when  machine-tuned  versions 
are  not  available. 


3.2  Preconditioners 


In  solving  linear  systems  Ax  = h using  iterative  techniques,  it  is  often  advantageous  to  precondition  the 
coefficient  matrix  A to  improve  convergence.  A preconditioner  M for  A needs  to  do  at  most  two  things: 
solve  the  system  Mx  = y , or  M'^x  ==  y.  Common  preconditioners  are  Jacobi  preconditioners,  where 
M = diag(A),  or  incomplete  factorizations  (LU,  or  Cholesky)  of  A. 

In  SparseLib-h-f  we  provide  preconditioner  class  structures  to  serve  this  purpose,  and  these  precondi- 
tioners can  be  used  with,  for  example,  the  IML-I--I-  Iterative  Methods  Library  (see  [5]).  The  current 
library  contains  the  following  preconditioners:  DiagPreconditionerO , for  diagonal  preconditioning, 
ICPreconditionerO,  for  incomplete  Cholesky  preconditioning,  and  ILUPreconditionerO,  for  incom- 
plete LU  preconditioning.  Each  preconditioner  provides  solveO  and  transpose-solveO  functionality, 
so  they  can  be  used  interchangeably  in  the  same  base  iterative  method  code  if  it  is  templated  for  a pre- 
conditioner. For  details  on  preconditioner  construction  and  use,  see  the  man  pages  in  appendix  A. 


4 File  I/O 


To  enable  the  use  of  SparseLib-f-f  with  previously  generated  sparse  matrices,  and  to  archive  newly 
created  sparse  matrices  for  later  use,  the  SparseLib-f- 1-  library  contains  several  functions  for  reading 
from  and  writing  to  plain  text  files  and  Harwell-Boeing  formatted  files.  (For  a detailed  description  of  the 
Harwell-Boeing  format,  see  [4].)  Note  that  although  the  internal  storage  scheme  indexing  in  SparseLib-|--f 
is  0-based,  file  I/O  is  performed  with  the  1-based  indexing  convention  to  provide  compatibility  with 
external  Fortran  libraries  and  other  packages  treating  sparse  matrices. 
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lunsyrametric  matrix  from  pores 


pores_l 


59 

2 

12 

45 

0 
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59 
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91 

101 

105 

115 

119 
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140 
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151 
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11 
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14 
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4 

5 

6 

13 

14 

3 

4 

5 

6 

7 

8 

15 

16 

5 

6 

7 

8 

15 

16 

5 

6 

7 

8 

9 

10 

17 

18 

7 

8 

10 

18 

7 

8 

9 

10 

19 

20 

9 

10 

19 

20 

1 

2 

11 

12 

13 

14 

21 

22 

11 

12 

13 

14 

21 

22 

3 

4 

11 

12 

13 

14 

15 

16 

23 

24 

13 

14 

16 

24 

5 

6 

13 

14 

15 

16 

17 

18 

25 

26 

15 

16 

18 

26 

7 

8 

15 

16 

17 

18 

19 

20 

27 

28 

17 

18 

20 

28 

9 

10 

17 

18 

19 

20 

29 

30 

19 

20 

29 

30 

11 

12 

21 

22 

23 

24 

21 

22 

24 

13 

14 

21 

22 

23 

24 

25 

26 

23 

24 

26 

15 

16 

23 

24 

25 

26 

27 

28 

25 

26 

28 

17 

18 

25 

26 

27 

28 

29 

30 

27 

28 

30 

19 

20 

27 

28 

29 

30 

29 

30 

-0.9481011349e+03  -0. 
0.9462545992e+03  0. 
-0.3005164596e+04  0. 
0.4731272996e+01  0. 
0.1552207555e+02  0. 


7178501646e+07 

7134130875e+07 

1293434629e+08 

3567021095e+05 

1628780334e+05 


0.4731272996e+01 

0.2334969309e+05 

-0.3680715121e+04 

-0.3120860678e+04 

0.3104415110e+04 


0.3574261854e+05 

-0.2461341087e+08 

0.6149543185e+07 

-0.3250082045e+07 

0.31914882106+07 


Figure  2:  A fragment  of  an  input  file  from  the  Harwell- Boeing  Sparse  Matrix  Collection. 


4.1  Reading  from  and  Writing  to  a Plain  Text  File 


CompCol_Mat_double  A; 

readtxtf ile_mat("input_file".  A) ; 

cout  » A; 

writetxtf ile_mat("output_f ile".  A) ; 


//or  CompRow_Mat_double , or  Coord_Mat_double 
//  read  from  file 
//  write  to  standard  out 
//  write  to  file 


Using  standard  output  (cout),  sparse  matrices  are  written  in  a plain  text  coordinate  format,  one  nonzero 
element  per  line,  with  each  line  containing  an  integer  row  index,  an  integer  column  index  and  a double 
value,  with  spaces  separating  the  three  fields.  To  ensure  that  the  row  and  column  dimensions,  M and 
N,  of  the  matrix  can  be  determined  implicitly  from  written  information,  we  include  in  the  output  the 
A(M-1,N-1)  matrix  element,  even  if  it  falls  outside  of  the  sparsity  pattern.  (This  convention  is  used  for 
sparse  matrix  saves  and  loads  in  MATLAB.)  The  functions  readtxtf  ile_mat  ( ) and  writetxtfilejnatO 
provide  file  access  for  reading  and  writing  sparse  matrix  information  in  this  plain  coordinate  format. 
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4.2  Reading  a Harwell-Boeing  Formatted  File 


void  readHB_iiif ©(const  chax  *filenaine,  int  &M,  int  &N,  int  &nz,  int  &nrh.s, 
int  verbose  = 0) 

void  readHB .header (FILE  *in_file,  char  *Title,  chax  *Key,  char  *mat_t3rpe, 

int  &Nrow,  int  ftNcol,  int  ftNnzero,  int  &Nrhs,  chax  *Ptrfmt, 
chax  *Indfmt,  chax  ♦Valfmt,  char  *Rhsfmt,  int  fePtrcrd, 
int  ftlndcrd,  int  ftValcrd,  int  ftRhscrd,  char  *Rhst3rpe) ; 


The  readHB  J.nfo  function  opens  and  reads  the  numerical  header  information  from  the  specified  Harwell- 
Boeing  file  and  returns  the  number  of  rows  and  columns  in  the  stored  matrix  (M  and  N),  the  number  of 
nonzeros  in  the  matrix  (nz),  and  the  number  of  right-hand-sides  stored  along  with  the  matrix  (nxhs). 
The  optional  verbose  parameter,  if  set  to  1,  sends  to  standard  output  more  detailed  information  from 
the  header,  including  title,  etc.  To  retrieve  and  save  to  variables  all  descriptive  header  information  from 
the  file  (such  as  “title”  or  “key’),  the  readHBJheader  is  available. 


void  readHB.mat (const  char  *fileneune, 
void  readHB.mat (const  chax  *filename, 
void  readHB_mat( const  char  ^filename. 


Coord_Mat_double  &A) 
CompCol_Mat_double  &A) 
CompRow_Mat_double  &A) 


The  readHB_mat  function  opens  and  reads  the  specified  file,  interpreting  its  contents  as  a sparse  matrix 
stored  in  the  Harwell/Boeing  standard  format  and  creating  a sparse  matrix  object  of  the  type  indicated 
in  the  calling  sequence. 


void  readHB.rhs (const  chax  *filenaine. 


MV_Vector_double  &b,  int  j=0) 


void  readHB.rhs (const  chax  *filename, 


MV.ColMat .double  &B) 


The  first  readHBxrhs  function  opens  and  reads  the  specified  file,  returning  a right-hand-side  vector  b.  If 
the  file  provides  a matrix  of  right-hand-sides,  (that  is,  a sequence  of  right-hand-side  vectors),  the  optional 
argument  j can  be  used  to  indicate  which  right-hand-side  is  desired.  The  default  reads  the  first  stored 
right-hand-side  (the  0th  column  of  the  right-hand-side  matrix).  The  second  form  of  the  readHBxrhs 
function  is  used  to  read  in  an  entire  multiple  right-hand-side  matrix,  assigning  it  to  B. 

4.3  Writing  a Sparse  Matrix  in  Harwell-Boeing  Format 


void  writeHB (const  chax  *filenaine,  const  Coord.Mat. double  A, 
const  chax  * Title,  const  char  * Key  ) 
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void  writeHB(const  chax  ^filename,  const  CompRovj_Mat_double  A, 
const  char  * Title,  const  char  * Key  ) 

void  writeHB( const  char  *filenajne,  const  CompCol_Mat_donble  A, 
const  char  * Title,  const  char  * Key  ) 

void  writeHB(const  char  *filenaine,  const  Coord_Mat_doTible  A) 

void  writeHB( const  char  ^filename,  const  CoinpRow_Mat_double  A) 

void  writeHB(const  char  *filenaine,  const  CompCol_Mat_donble  A) 


The  writeHB  function  opens  the  named  file  and  writes  the  specified  matrix  to  that  file  in  Harwell- 
Boeing  format.  If  the  Title  and  Key  arguments  are  not  supplied,  default  values  are  used  to  indicate 
that  SparseLib++  generated  the  file. 


5 Programming  Considerations 

5.1  Integrating  with  Fortran 

Since  SparseLib++  matrices,  like  C and  C++  arrays,  have  0-based  indexing,  some  care  must  be  used 
to  ensure  compatibility  with  Fortran  subroutines.  Note  that  this  is  an  issue  only  with  sparse  matrices, 
since  the  specific  row/column  index  values  are  part  of  the  data. 

There  are  several  ways  to  handle  this  issue:  (1)  one  could  create  a copy  of  the  matrix  with  each  row 
and/or  column  index  incremented  by  one,  (2)  make  SparseLib++  matrices  explicitly  1-based  (or  use 
arbitrary  bases),  or  (3)  modify  the  0-based  arguments  to  1-based  subroutines. 

Solution  (1)  is  easy  to  implement  but  rather  expensive  in  practice,  requiring  an  extra  copy  of  the  matrix 
as  well  as  a 0{nz)  algorithm  to  update  the  indices.  Solution  (2)  seems  reasonable,  particularly  for 
Fortran  programmers,  but  makes  SparseLib++  incompatible  with  native  C arrays.  There  are  several 
other  inconsistencies  brought  up  by  violating  the  basic  indexing  scheme  of  C.  Solution  (3)  can  sometimes 
be  used,  but  is  not  a universal  solution.  For  example,  the  Fortran  Sparse  BLAS  matrix-vector  multiply 
routine  (y  <—  aAx-\-fiy)  for  a Coordinate  storage  sparse  matrix  and  0-based  C/C++  array  can  be  called 
as: 


DC00MH_("N",  M,  N,  1,  alpha.  A,  &x(l),  Idb,  beta,  &y(l),  N, 
work,  iwork) ; 

Similarly,  for  a more  general  matrix-matrix  multiply  , we  can  write  (C  <—  aAB  + /3C) 
DC00HM_(’'N",  M,  K,  K,  alpha.  A,  &B(i,j)+  1,  Idb,  beta,  &C(k,l)  + l,  Idc, 
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work,  iwork) ; 


This  scheme  makes  assumptions  about  how  the  0-based  parameters  are  used  internally,  and  therefore  its 
correctness  cannot  be  guaranteed. 

Note,  also,  that  an  expression  such  as  "&x(l)"  assumes  a vector  has  at  least  two  elements.  If  this 
causes  an  error  when  bounds  checking  is  turned  on,  the  equivalent  expression  "&xCO]+l"  can  be  used  to 
circumvent  the  bounds  checking  in  this  case. 


5.2  Ordering  of  Sparse  Matrix  Elements 

Many  of  the  sparse  matrix  formats  do  not  require  a unique  ordering  of  their  elements.  Most  obvious 
in  this  category  is  coordinate  storage,  but  other  (more  structured)  formats  share  this  ambiguity.  The 
compressed  row  storage  format,  for  example,  cissumes  no  explicit  ordering  among  the  elements  in  each 
row,  but  there  are  varying  levels  of  ordering  which  may  be  present,  and  advantages  from  each  level.  No 
assumed  ordering  is  the  most  general;  ordering  diagonal  elements  from  each  row  first  within  each  row  is 
convenient  for  scaling  and  Jacobi  preconditioners;  ordering  within  rows  by  column  index  is  probably  the 
most  natural  and  provides  for  most  efficient  addressing.  Although  the  current  release  of  SparseLib++ 
does  not  assume  ordering,  future  releases  will  be  equipped  to  order  the  elements  to  suit  user  requirements. 
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Coord_Mat_double 


Name  Coord-Mat.double 

Declaration  #include  "coord. h" 

class  CoordJVIat -double 


Description  A coordinate  formatted  sparse  matrix  class.  Access  to  matrix  elements  as  A(i,  j) 
are  0-based  (i.e.  A(0,0)  is  the  first  element). 

The  nonzero  matrix  elements  and  their  indices  are  stored  in  three  vectors,  with 
val  ( ) holding  the  element  values,  row_ind( ) holding  the  row  indices,  and  col_iiid( ) 
holding  the  column  indices.  For  example,  a Coord_Mat_double  object  specifying 


the  matrix 

/ 

1 

2 

0 

0 

3 

\ 

4 

5 

6 

0 

0 

0 

7 

8 

0 

9 

0 

0 

0 

10 

0 

V 

11 

0 

0 

0 

12 

/ 

may  have  the  following  vectors 

in 

internal 

storage: 

valO 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

rewind  0 

0 

0 

0 

1 

1 

1 

2 

2 

2 

3 

4 

4 

col-indO 

0 

1 

4 

0 

1 

2 

1 

2 

4 

3 

0 

4 

Constructors /Destructors 

Coord-Mat-double  ( ) 

Construct  a null  0x0  matrix. 

Coord-Mat-double  ( const  Coord_Mat_double  SzC  ) 

Construct  a copy  of  the  coordinate  matrix  C. 

Coord-Mat-double  ( const  CompColJVIat-double  icC  ) 

Construct  a coordinate  matrix  from  a given  compressed  column  representation. 

Coord-Mat-double  ( const  CompRow-Mat_double  &C  ) 

Construct  a coordinate  matrix  from  a given  compressed  row  representation. 

Coord-Mat-double  ( int  M,  int  N,  int  nz,  double  *val,  int  *r,  int  *c  ) 

Construct  a coordinate  matrix  of  size  M x N,  with  nz  nonzeros,  using  the  values 
given  in  val  []  , and  row  and  column  indices  given  in  r []  and  c [] . The  vectors  r 
and  c are  assumed  to  be  0— based. 
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~CoordJVIat_double  ( ) 

Matrix  destructor. 

Access  and  Information 

double  operator  ( int  i,  int  j ) 

Return  A{i,j).  (Returns  zero  if  matrix  element  not  found  in  sparse  structure.) 
doubled  set  ( int  i,  int  j ) 

Assign  A{i,j)  a value.  Reports  an  error  to  stderr  and  exits  the  program  if  the 
assignment  violates  the  sparsity  structure  of  A (i.e.  causes  fill-in).  For  dynamically 
growing  a sparse  matrix,  use  an  appropriate  class. 

doubled  val  ( int  i ) 

Return  the  ith  element  of  the  nonzero  value  storage  vector, 
int  row_ind  ( int  i ) 

Return  the  row  index  of  the  element  stored  in  val(i). 
int  colJnd  ( int  i ) 

Return  the  column  index  of  the  element  stored  in  val(i). 
int  dim  ( int  i ) 

Return  the  size  of  the  matrix  along  dimension  i. 
int  size  ( int  i ) 

Return  the  size  of  the  matrix  along  dimension  i.  (Same  as  dim().) 
int  NumNonzeros  ( ) 

Return  the  number  of  nonzeros  in  the  matrix. 

Standard  Output 

friend  ostream&:  operator<<  ( ostream  &os,  const  Coord_Mat_double  &C  ) 
Print  nonzero  matrix  elements  one  per  line  in  the  format  <i>  < j > < val  > . 

See  also  CompRow_Mat_double,  CompCol_Mat_double,  MV.Vector  (MV-|— 1-), 
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Name  CompRow_Mat_double 

Declaration  #include  "comprow.h" 

class  CompRowJVIat -double 

Description  A compressed  row  formatted  sparse  matrix  class.  Access  to  matrix  elements  as 
A(i,j)  are  0-based  (i.e.  A(0,0)  is  the  first  element). 

The  nonzero  elements  and  index  information  are  stored  in  three  vectors,  with  valO 
holding  the  values  of  the  nonzero  elements,  row_ptr()  holding  pointers  to  the  first 
element  in  each  row,  and  col-i.nd( ) holding  a column  index  for  each  of  the  elements 
in  val().  An  additional  element  is  appended  to  the  row_ptr()  array  specifying  the 


number  of  nonzero  array  elements. 

For 

example,  a CompRow -Mat-double  object 

specifying  the  matrix 

/ 1 

2 

0 

0 

3 ^ 

4 

5 

6 

0 

0 

0 

7 

8 

0 

9 

0 

0 

0 

10 

0 

V 11 

0 

0 

0 

12  yi 

may  have  the  following  vectors  in  internal  storage: 


row-ptr  ( ) 

0 

3 

6 

9 

10 

12 

val() 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

coU.nd() 

0 

1 

4 

0 

1 

2 

1 

2 

4 

3 

0 

4 

Constructors /Destructors 

CompRow -Mat-double  ( ) 

Construct  a null  0x0  matrix. 

CompRo-w -Mat-double  ( const  CompRow -Mat-double  &R  ) 

Create  a copy  of  the  compressed  row  matrix  R. 

CompRowr -Mat-double  ( const  CompCoord_Mat_double  &R  ) 

Construct  a compressed  row  matrix  from  a given  coordinate  representation. 

CompRo-w -Mat-double  ( const  CompCol_Mat-double  &R  ) 

Construct  a compressed  row  matrix  from  a given  compressed  column  representation. 

CompRow -Mat -double  ( int  M,  int  N,  int  nz,  double  *val,  int  *r,  int  *c  ) 

Construct  a compressed  row  matrix  of  size  M x N,  with  nz  nonzeros,  using  the  values 
given  in  val  [] , and  row  pointers  and  column  indices  given  in  r []  and  c [] . The 
vectors  r and  c are  assumed  to  be  0— based. 
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~CompRow -Mat -double  ( ) 

Matrix  destructor. 

Access  and  Information 

double  operator  ( int  i,  int  j ) 

Return  A{i,j).  (Returns  zero  if  matrix  element  not  found  in  sparse  structure.) 
doubled  set  ( int  i,  int  j ) 

Assign  A{i,j)  a value.  Reports  an  error  to  stderr  and  exits  the  program  if  the 
assignment  violates  the  sparsity  structure  of  A (i.e.  causes  fill-in).  For  dynamically 
growing  a sparse  matrix,  use  an  appropriate  class. 

doubled  val  ( int  i ) 

Return  the  ith  element  of  the  nonzero  value  storage  vector, 
int  row-ptr  ( int  i ) 

Return  the  row  pointer  associated  with  row  i. 
int  col-ind  ( int  j ) 

Return  the  column  index  of  the  element  stored  in  val(j). 
int  dim  ( int  i ) 

Return  the  size  of  the  matrix  along  dimension  i. 
int  size  ( int  i ) 

Return  the  size  of  the  matrix  along  dimension  i.  (Same  as  dimO.) 
int  NumNonzeros  ( ) 

return  the  number  of  nonzeros  in  the  matrix. 

Standard  Output 

friend  ostream&:  operator<<  ( ostream  &:os, 

const  CompRow -Mat-double  &:R  ) 

Print  matrix  elements  one  per  line  in  the  format  <i>  < j > < val  > . 

See  also  CompCol-Mat.double,  CoordJMat-double,  MV.Vector  (MV-t— h), 

SparseLib-f-f-  File  I/O 
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Name  CompCol_Mat-double 

Declaration  #include  "compcol.h" 

clciss  CompCol_Mat_double 

Description  A compressed  column  formatted  sparse  matrix  class.  Access  to  matrix  elements  as 
A(i,  j)  are  0-based  (i.e.  A(0,0)  is  the  first  element). 

The  nonzero  elements  and  index  information  are  stored  in  three  vectors:  val() 
holding  the  values  of  the  nonzero  elements,  col_ptr()  holding  pointers  to  the  first 
element  in  each  column,  and  row_ind()  holding  a row  index  for  each  of  the  elements 
in  val().  An  additional  element  is  appended  to  the  col_ptr()  array  specifying  the 
number  of  nonzero  array  elements.  For  example,  a CompColJVIat.double  object 
specifying  the  matrix 

/ 1 2 0 0 3 \ 

4 5 6 0 0 

0 7 8 0 9 

0 0 0 10  0 

\ 11  0 0 0 12  / 

may  have  the  following  vectors  in  internal  storage: 


col-ptr() 

0 

3 

6 

8 

9 

12 

valO 

1 

4 

11 

2 

5 

7 

6 

8 

10 

3 

9 

12 

row-indO 

0 

1 

4 

0 

1 

2 

1 

2 

3 

0 

2 

4 

Constructors /Destructors 

CompCol-Mat-double  ( ) 

Construct  a null  0x0  matrix. 

CompCol_Mat_double  ( const  CompCol-Mat.double  &;C  ) 

Create  a copy  of  the  compressed  column  matrix  C. 

CompColJMat -double  ( const  CompCoord_Mat_double  &C  ) 

Construct  a compressed  column  matrix  from  a given  coordinate  representation. 

CompColJMat-double  ( const  CompRow-Mat_double  &:C  ) 

Construct  a compressed  column  matrix  from  a given  compressed  row  representation. 

CompColJMat-double  ( int  M,  int  N,  int  nz,  double  *val,  int  *r,  int  *c  ) 

Construct  a compressed  column  matrix  of  size  M x N,  with  nz  nonzeros,  using  the 
values  given  in  val  [] , and  row  indices  and  column  pointers  given  in  r []  and  c [] . 
The  vectors  r and  c are  assumed  to  be  0— based. 
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~CompColJVIat_doubie  ( ) 

Matrix  destructor. 

Access  and  Information 

double  operator  ( int  i,  int  j ) 

Return  A{i,j).  (Returns  zero  if  matrix  element  not  found  in  sparse  structure.) 
doubled  set  ( int  i,  int  j ) 

Assign  A{i,  j)  a value.  Reports  an  error  to  stderr  and  exits  the  program  if  the 
assignment  violates  the  sparsity  structure  of  A (i.e.  causes  fill-in).  For  dynamically 
growing  a sparse  matrix,  use  an  appropriate  class. 

doubled  val  ( int  z ) 

Return  the  ith  element  of  the  nonzero  value  storage  vector, 
int  row  Jnd  ( int  i ) 

Return  the  row  index  of  the  element  stored  in  val(i). 
int  col-ptr  ( int  j ) 

Return  the  column  pointer  associated  with  column  j. 
int  dim  ( int  i ) 

Return  the  size  of  the  matrix  along  dimension  i. 
int  size  ( int  i ) 

Return  the  size  of  the  matrix  along  dimension  i.  (Same  as  dim().) 
int  NumNonzeros  ( ) 

return  the  number  of  nonzeros  in  the  matrix. 

Standard  Output 

friend  ostream&  operator<<  ( ostream  &;os, 

const  CompCol-Mat-double  izC  ) 

Print  nonzero  matrix  elements  to  standard  output,  one  per  line,  in  the  format  <i> 
< j ><  val  > . 

See  also  CompRow_Mat_double,  Coord_Mat_double,  MV.Vector  (MV-I— 1-), 

SparseLib-|--f  File  I/O 


Version  1.5 


15 


June  27,  1996 


DiagPreconditioner 


SparseLib++ 


DiagPreconditioner 


Name 

DiagPreconditioner  — Diagonal  Preconditioner  Class 

Declaration 

^include  "diagpre.h" 

class  DiagPreconditioner 

Description 

Implements  diagonal  preconditioning  for  use  with  IML++  iterative  methods.  Presently, 
implementations  for  SparseLib++  compressed  row  and  compressed  column  matrix 
formats  {CompColJiJai-double  and  Comp  Row -Mat-iouhle,  respectively)  have  been 
provided.  In  general,  users  of  IML-f  + will  not  need  to  access  member  functions  of 
this  class  except  to  create  an  instance  of  a preconditioner. 

Constructors /Destructors 


DiagPreconditioner  ( const  CompCol_Mat_double&  A ) 

Construct  a diagonal  preconditioner  from  the  matrix  A. 

DiagPreconditioner  ( const  CompRowJMat_double&:  A ) 

Construct  a diagonal  preconditioner  from  the  matrix  A. 

~DiagPreconditioner  ( void  ) 

Reclaim  memory  space. 

Member 

Functions 

Vector-double  solve  ( const  Vector_double&:  h )const 

Perform  the  preconditioning,  that  is,  return  the  solution  of  the  linear  system  with 
the  preconditioner  and  the  vector  h. 

Vector-double  trans-solve  ( const  Vector_double&:  h )const 

Perform  the  transpose  preconditioning,  that  is,  return  the  solution  of  the  linear 
system  with  the  transposed  preconditioner  and  the  vector  h.  For  the  diagonal  pre- 
conditioner  (which  is  trivially  symmetric),  this  is  the  same  as  the  solve()  member 
function. 

Example 

For  examples  of  the  use  of  this  class,  see  the  examples  provided  with  the  descriptions 
of  the  IML-|--(-  iterative  method  functions. 

See  Also 

SparseLib-f-l- 

Version  1.5 


16 


June  27,  1996 


ICPreconditioner 


SparseLib++ 


ICPreconditioner 


Name 

ICPreconditioner  — Incomplete  Cholesky  Preconditioner 

Declaration 

#include  "icpre.h" 

class  ICPreconditioner 

Description 

Implements  incomplete  Cholesky  preconditioning  for  use  with  IML++  iterative 
methods.  Presently,  implementations  for  SparseLib++  compressed  column  and 
compressed  row  matrix  formats  [Comp Col-Mat. double  and  CompRow-Mai.double, 
respectively)  have  been  provided.  In  general,  users  of  IML++  will  not  need  to  access 
member  functions  of  this  class  except  to  create  an  instance  of  a preconditioner.  The 
matrix  A must  be  symmetric  positive  definite. 

The  present  implementation  of  incomplete  Cholesky  factorization  does  not  create 
any  fill  for  structural  zero  elements  (i.e.,  it  is  an  implementation  of  IC(0)),  but  it 
does  modify  non-zero  elements. 

Constructors /Destructors 


ICPreconditioner  ( const  CompCol_Mat_double&  A ) 

Construct  an  incomplete  Cholesky  preconditioner  from  the  matrix  A. 

ICPreconditioner  ( const  CompRowJMat_double&;  A ) 

Construct  an  incomplete  Cholesky  preconditioner  from  the  matrix  A. 

~ICPreconditioner  ( void  ) 

Reclaim  memory  space. 

Member 

Functions 

Vector-double  solve  ( const  Vector_double&  b )const 

Perform  the  preconditioning,  that  is,  return  the  solution  of  the  linear  system  with 
the  preconditioner  and  the  vector  b. 

Vector-double  trans_solve  ( const  Vector-double&;  b )const 

Perform  the  transpose  preconditioning,  that  is,  return  the  solution  of  the  linear  sys- 
tem with  the  transposed  preconditioner  and  the  vector  b.  For  the  IC  preconditioner 
(which  is  symmetric),  this  is  the  same  as  the  solve()  member  function. 

Example 

The  following  example  program  uses  IML-t— (-  in  conjunction  with  SparseLib-|--|-  to 
solve  a linear  system  with  CG.  The  program  reads  in  a matrix  and  right-hand  side 
stored  in  Harwell-Boeing  format  from  the  file  test.dai.  An  initial  guess  of  0 is  made 
for  the  solution  and  the  system  is  solved  using  CG  and  an  incomplete  Cholesky 
preconditioner. 
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#iiiclude  <iostreain.h> 
#include  <stdlib.h> 
#include  "cg.h" 
#include  "icpre.h" 
#include  "compcoll .h" 
#include  "iohb.h" 
#include  "vector.h" 
tfinclude  "blasl.h" 


int  mainO 

double  tol  = l.e-6; 
int  result,  maxit  = 150; 


CompCol_Mat_double  A; 
readHB ("test .hb".  A); 

Vector_double  x(A.dim(l),  0.0),  b; 
readHB ("test .hb",  b) ; 

ICPreconditioner  D(A) ; 

result  = CG(A,  x,  b,  D,  maxit. 


//  Create  a matrix 
//  Read  matrix  data 

//  Solution,  rhs  vectors 
//  Read  rhs  data 

preconditioner 

//  Solve  system  with  CG 


//  IC 
tol)  ; 


cout  « "CG  flag  = " <<  result  « endl; 
cout  <<  "iterations  performed:  " <<  maxit  <<  endl; 
cout  <<  "tolerance  achieved  : " <<  tol  <<  endl; 
cout  <<  x; 


return  0; 

} 


See  Also  CG 

SparseLib-|-+ 

R.  Barrett  ET  al.,  Templates  for  the  Solution  of  Linear  Systems:  Building  Blocks 
for  Iterative  Methods,  SIAM  Press,  Philadelphia,  1994. 

J.  Meijerink  and  H.  a.  van  DER  Vorst,  An  iterative  solution  method  for  linear 
systems  of  which  the  coefficient  matrix  is  a symmetric  M -matrix,  Math.  Comp., 
31  (1977),  pp.  148-162. 
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ILUPreconditioner 


Name 

ILUPreconditioner  — Incomplete  LU  Preconditioner 

Declaration 

#include  "ilupre.h" 

class  CompCoULUPreconditioner 
class  CompRow -ILUPreconditioner 

Description 

Implement  incomplete  LU  preconditioning  for  use  with  IML++  iterative  meth- 
ods. Presently,  implementations  for  SparseLib-|— f compressed  row  and  compressed 
column  matrix  formats  [CompCol-Mai.double  and  CompRow-Mat-douhle,  respec- 
tively) have  been  provided.  In  general,  users  of  IML-f-l-  will  not  need  to  access 
member  functions  of  this  class  except  to  create  an  instance  of  a preconditioner. 

This  implementation  of  incomplete  LU  factorization  does  not  create  any  fill  for 
structural  zero  elements  (i.e.,  it  is  an  implementation  of  ILU(O)),  but  it  does  modify 
non-zero  elements. 

Constructors/Destructors 


CompCoULUPreconditioner  ( const  CompColJV[at_double&:  A ) 

Construct  an  incomplete  LU  preconditioner  from  the  matrix  A. 

~CompCoULUPreconditioner  ( void  ) 

Reclaim  memory  space. 

CompRow -ILUPreconditioner  ( const  CompRowJVlat_double&  A ) 

Construct  an  incomplete  LU  preconditioner  from  the  matrix  A. 

~CompRowJ[LUPreconditioner  ( void  ) 

Reclaim  memory  space. 

Member 

Functions 

Vector-double  solve  ( const  Vector-double&  b )const 

Perform  the  preconditioning,  that  is,  return  the  solution  of  the  linear  system  with 
the  preconditioner  and  the  vector  b. 

Vector-double  trans_solve  ( const  Vector_double&  b )const 

Perform  the  transpose  preconditioning,  that  is,  return  the  solution  of  the  linear 
system  with  the  transposed  preconditioner  and  the  vector  b. 

Example 

The  following  example  program  uses  IML-t— t-  in  conjunction  with  SparseLib-|— f to 
solve  a linear  system  with  GMRES.  The  program  reads  in  a matrix  and  right-hand 
side  stored  in  Harwell- Boeing  format  from  the  file  test.dat.  An  initial  guess  of  0 is 
made  for  the  solution  and  the  system  is  solved  using  GMRES  and  an  incomplete 
LU  preconditioner. 
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ttinclude 
#include 
#iiiclude 
#include 
# include 
# include 
#include 
#include 


<iostream.h> 
<stdlib.h> 
"gmres .h" 
"ilupre .h" 
"compcoll .h" 
"iohb.h" 
"vector .h" 
"blasl.h" 


int  mainO 

{ 

double  tol  = l.e-6; 

int  result,  m = 32,  maxit  = 150; 

CompCol_Mat .double  A; 
readHBC'test  .hb".  A); 

Matrix.double  H(m+1,  m,  0.0); 
Vector.double  x(A.dim(l),  0.0),  b; 
readHBC'test  .hb",  b)  ; 

CompCol.ILUPreconditioner  D(A) ; 

result  = GMRES (A,  x,  b,  D,  H,  m 


//  Create  a matrix 
//  Read  matrix  data 

//  H matrix 

//  Solution,  rhs  vectors 
//  Read  rhs  data 

//  ILU  preconditioner 

maxit,  tol);  //  Solve  system  with  GMRES 


cout  <<  "GMRES  flag  = " <<  result  <<  endl; 
cout  « "iterations  performed:  " « maxit  « endl; 
cout  <<  "tolerance  achieved  : " <<  tol  <<  endl; 
cout  « x; 


return  0; 

} 


See  Also  GMRES 

SparseLib++ 

R.  Barrett  ET  al.,  Templates  for  the  Solution  of  Linear  Systems:  Building  Blocks 
for  Iterative  Methods,  SIAM  Press,  Philadelphia,  1994. 

J.  Meijerink  and  H.  a.  van  DER  Vorst,  An  iterative  solution  method  for  linear 
systems  of  which  the  coefficient  matrix  is  a symmetric  M -matrix.  Math.  Comp., 
31  (1977),  pp.  148-162. 
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